Spontaneous formation of synchronization clusters in homogenous neuronal ensembles 

induced by noise and interaction delays 
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Spontaneous formation of clusters of synchronized spiking in a structureless ensemble of equal 
stochastically perturbed excitable neurons with delayed coupling is demonstrated for the first time. 
The effect is a consequence of a subtle interplay between interaction delays, noise and the excitable 
character of a single neuron. Dependence of the cluster properties on the time-lag, noise intensity 
and the synaptic strength is investigated. 



Collective behavior in large ensembles of physiological 
and inorganic systems can be reduced to that of coupled 
oscillators engaged in various synchronization phenom- 
ena. In terms of macroscopic coherent rhythms, it may 
either be the case where all the units are recruited into a 
giant component or the case of cluster states character- 
ized by exact or in-phase intra-subset and lag inter-subset 
synchronization. The spontaneous onset of cluster states 
is of particular interest to neuroscience [l[ for the con- 
jectured role in information encoding, as well as for par- 
ticipating in motor coordination or accompanying some 
neurological disorders. The approach to clustering has 
mostly relied on modeling neurons as autonomous oscil- 
lators, treating separately the question of whether the 
proposed mechanisms may be robust under noise Q and 
transmission delays [3|- We explore a new mechanism 
which rests on the excitable character of neuronal dy- 
namics and mutual adjustment between noise and time 
delay to yield the self-organization into functional mod- 
ules within an otherwise unstructured network. 

For the instantaneous couplings, the research on pop- 
ulations of excitable neurons has covered pattern forma- 
tion due to local inhomogeneity [4] , or has invoked a sce- 
nario where noise enacts a control parameter tuning the 
dynamics of ensemble averages between the three generic 
global regimes [5] . Distinct from the layout with complex 
connection topologies, here it is demonstrated how cou- 
pling delays do alter the latter landscape in a significant 
fashion, giving rise to an effect one may dub the clus- 
ter forming time-delay- induced coherence resonance. In 
part, the strategy to analyze global dynamics rests on de- 
riving the mean-field (MF) approximation for the exact 
system. The likely gain from the MF treatment is at least 
twofold: except for allowing one to extrapolate what oc- 
curs in the thermodynamic limit N — » oo, it may serve as 
an auxiliary means to discriminate between the effects of 
noise and time delay. Unexpectedly, the MF model un- 
dergoes a global bifurcation at certain parameter values 



where the exact system shows an onset of cluster states. 

Network dynamics and the tools to analyze it- We fo- 
cus on an TV-size population of all-to-all diffusively cou- 
pled Fitzhugh-Nagumo neurons, whose local dynamics is 
set by 

N 

edxi = (xi - x 3 j3 - y + I)dt + ^[xj(t - r) - Xi(t)]dt, 



dyi = (x + b)dt + V2DdWi, 



(1) 



where the activator variables xi embody the membrane 
potentials, while the recovery variables yi mimic the ac- 
tion of the K + membrane gating channels, c denotes 
the synaptic strength and r stands for the coupling de- 
lay, both parameters for simplicity assumed homoge- 
neous across the ensemble. The V2DdWi terms rep- 
resent stochastic increments of the independent Wiener 
processes, i.e. the white noise. For the external stimula- 
tion holds I = 0, whereas the small parameter e = 0.01 
warrants a clear separation between the fast and slow 
time scales. Selecting b = 1.05, the neurons are poised 
near the Hopf bifurcation threshold 6 = 1, which places 
them in an excitable regime where each possesses a sin- 
gle equilibrium. An adequate stimulation, be it by the 
noise or the interaction term, may evoke a large excur- 
sion of membrane potential, passing through the spiking 
and refractory states before it loops back to rest. 

To characterize the degree of correlation be- 
tween the firing events, we use primarily the 
interneuron spike train coherence [6] = 
Y%=iMk)X 3 (k)/ ^Y,lUMk)X 3 (k). This requires 
one to split the simulation period T into bins k of 
length A = T/m, awarding each neuron a variable 
Xi(k) = 1(0), contingent on whether a spike was 
triggered or not within the given time bin, respectively. 
As with all the quantities below, we have been careful 
to exclude from calculations the transient behavior. The 
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spike threshold and the time bin are set to Xq = 1 
and A = 0.008, verifying that no change of the results 
occurred if Xq or A were reduced. The distribution 
of the Kij values may serve to distinguish between the 
homogeneous and clustered network states. Another 
aspect we are interested in is whether the clustered 
states are monostable or coexistent with the homoge- 
neous ones at the given network size. To probe this, 
we have monitored if the values of the global coherence 

N 

K = N(N—i) 5^ K ij f° r different realizations at the 

fixed parameters clumped together, expecting bunching 
into distinct groups as evidence of mult ist able behavior. 

Addressing the temporal structure of the network 
states, it is useful to look into the distribution of the 
local neuron jitters t\ [7]. They represent the normal- 
ized variations of the interspike intervals extracted 
from Xi (t), n = v / < T k>~< T k >7 < T k >, with 
smaller values indicating better regularity. The modality 
and the width of the Ti distribution over the population 
may serve as rough indices on how the cluster dynamics 
is mutually adjusted. In the final part, we analyze the 
behavior of the ensemble averages X = 1/N Y^iLi x i an d 
Y = l/N^2 i=1 yi, where the former increases if a larger 
fraction of neurons fire in synchrony. The results for the 
exact system are compared to those of the approximate 
MF model [9[ . The latter presents a two-dimensional set 
of delayed differential equations 



dt 



dY(t) 
dt 



M |i _ c _ X (t) 2 + y/[c - 1 + X(t) 2 ] 2 + 40} 
Y(t)+c[X(t-T)-X(t)}, 

: X(t) + b, (2) 



derived within a cumulant approach by employing the 
Gaussian approximation. 

We note that the results for the exact system refer to 
a network of N = 200 neurons, applying independently 
a method from [8] to verify no qualitative changes in the 
clustering behavior for larger N. 

Results- To get a sense of what may be the param- 
eter ranges to admit the cluster states, we plot the c- 
families of the k curves in dependence of D for different 
r. Without the delay, the curves would conform to a 
stereotype profile, where one distinguishes between the 
three "regular" segments for very small, intermediate and 
large D, showing first a reduced k due to incoherent os- 
cillations, then steady high values for the coherent ones 
and the decaying segment at D where the stochastic dy- 
namics prevails. However, from Fig. [T] we learn how 
this is upheld for some r, say r = 11, but is violated 
manifestly at the "cluster-resonant" values r = 2,6,10. 
The "wells" seen at approximately D G (0.001,0.003) 
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FIG. 1. (color online). Profiles of the k(D) families of curves 
over the synaptic strengths c = 0.08,0.1 and 0.12 display 
strong dependence on the delay, increasing from r = 2 in (a), 
r = 6 in (b), r = 10 in (c) to r = 11 in (d). The location of 
"wells" may point to the emergence of the clustered states. 

in Figs. QJb) and QJc) may occur for just two reasons, 
as k decreases either for the incoherent or the clustered 
states. The latter alternative is supported by the co- 
herence matrices in Fig. [3j which are discussed shortly. 
The importance of the D - r adjustment for the cluster- 
ing effect is also witnessed by the c-dependence within 
the families in Fig. [1] the stronger the interaction term, 
the more salient is the picture of "irregularity" sections 
immersed into a "regular" curve profile. Increasing the 
delay, the cluster states first occur, apparently monos- 
table, around r = 2 for the small D = 0.00025, whereby 
the typical phase portrait (PP) projection shows twisted 
orbits with two clearly discernable segments, see Fig. 
[2fa). These reflect the two macroscopic fractions of 
the population firing alternately, such that the homo- 
geneous network dynamically splits into clusters of mu- 
tually synchronized neurons, with the clusters locked in 
antiphase. The frequency entrainment is indicated by the 
shape of the T{ distribution, which peaks sharply around 
< t >m= 0.01. We tested the invar iance of clustering 
with N via the asymptotical behavior of the quantity 
Xn = Ty^-r, where a\ =< X(t) 2 > t - < X(t) >f 

and G 2 Xi =< Xi(t) 2 > t — < Xi(t) > 2 holds. If the clus- 
ter states endure, there should be a residual component 
x(oo) G (0,1) in the large N limit [8]. For this and 
the remaining cases, the onset of such a regime is found 
around N w 200, implying that no qualitatively novel 
phenomena occur above this system size. An interesting 
observation is that the cluster configuration {A/i,^}, 
determined by the fractions' sizes, fluctuates around the 
ratio 2 : 1 for different stochastic realizations and ap- 
pears to aggregate with enhancing N. For certain r, the 
2-cluster state also emerges outside the D-region delim- 
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FIG. 2. Global PPs for the 2-cluster states show twisted LCs, 
whereby the two discernable segments reflect the alternate 
firing of the neuron subsets. The N1/N2 ratio depends on the 
interplay of D and r, as seen from the examples r — 2, D — 
0.00025, c = 0.1 in (a) and r = 5, D = 0.0005, c = 0.1 in (b). 

iting the incoherent and coherent global regimes. This 
holds for r = 5 and D e (0.0004, 0.0008), where the clus- 
ter layout is also such that if one is active, the other re- 
mains refractory. The distribution maintains a narrow 
form, but its maximum shifts to<r> m ^0.19. Though 
one retrieves the general picture from above, a variance is 
that larger r seems to favor the partition N1/N2 ~ 1 : 1, 
see PP in Fig. [2fb). The 1 : 1 ratio is preferred both for 
increasing N and if the delay is set to r = 6. 

The clustered states so far may be cast as stationary in 
the sense of stability against neurons switching between 
the clusters. We also report on the existence of 3-cluster 
states that may be considered " dynamical" , with the neu- 
rons able to jump to and from clusters. Such an outcome 
arises for the stronger noise D « 0.0013, once the delay is 
increased to r = 10. To underline the difference between 
the stationary and dynamical clustered states at r = 5 
and r = 10, we plot side-by-side the corresponding pair- 
wise coherence matrices {/%}, see Figs. [3ja) and[3fb), 
where the network nodes have been rearranged by a hier- 
archical clustering algorithm according to a form of met- 
ric distance that has the most coherent nodes the closest. 
This makes it explicit how the intercluster coherence for 
the 2-cluster state is virtually negligible with respect to 
the 3-cluster case. Loosely speaking, within an unstable 
three-part population division, when a certain fraction 
is firing, the other is refractory and the neurons in the 
smallest cluster are at rest (excitable). This less clear 
separation is also apparent when comparing the nodal 
degree distributions in cases r = 5 and r = 10, obtained 
if one assumes {&ij} to provide weights for the network 
whose links stand for the correlated dynamics between 
the neurons. For r = 5, the bimodal degree distribution 
is clearly seen without raising the connectivity thresh- 
old, whereas for r = 10 the initially smeared three-modal 
distribution refines after some thresholding is performed. 
The rationale of dynamical clustering may best be un- 
derstood by analyzing the n distribution in the 3-cluster 
state. Apart from being wider than in the 2-cluster state, 
it peaks at a much smaller value < r > m ~ 0.09, im- 
plying the more regular neuron firing. For this to hold, 
synchrony within the clusters has to be of intermittent 
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FIG. 3. (color online). Rearranged coherence matrices for 
r = 5, D = 0.0005, c 0.1 in (a) and r = 10, D = 0.0013, c 
0.1 in (b) imply the strong cluster separation in the 2-cluster 
states and mixing between the clusters in the 3-cluster case. 
Darker shading reflects higher coherence, 
nature, such that the neurons once engaged in synchro- 
nized spiking are much more likely to do so again. 

Aa understanding of clustering mechanism is revealed 
by comparing the typical PPs of neurons participant in 
the homogeneous coherent state and the clustered states, 
see Figs. HJa) and|Hb). A striking feature in the latter 
case is a kink at the refractory branch of the slow mani- 
fold. The appearance of a kink is the key manifestation 
of the D — r co-effect, that consists in separating the 
ensemble into clusters and maintaining the proper phase 
difference between them. The purpose of the kink is to 
keep the neurons frustrated long enough at the refractory 
branch before being allowed to slide down to its left knee. 
This may be imagined as a form of a lock-and-release be- 
havior, where the delay primarily gives rise to the first, 
and noise to the second part. If a fraction of the en- 
semble were to move beyond the left knee and the other 
were to lag behind, the split should be amplified with 
each population cycle, eventually becoming resilient to 
perturbation precisely due to trapping at the refractory 
branch. For trapping to be successful, the kink has to 
be placed properly, approximately where the dynamics 
of the representative point is most susceptible to pertur- 
bation along the slow manifold. Then, for a brief pe- 
riod, due to an influence from a^, the evolution of yi is 
locally accelerated, becoming comparable to a speed of 
change in the direction orthogonal to the slow manifold, 
driven by the spiking fraction of the population. Note 
that the trapping interval has to be adjusted so that the 
entire population is entrained to a single frequency of fir- 
ing. The latter matches the one in delay-free case, which 
warrants stability against perturbations. The arguments 
above and the numerical data seem to indicate how the 
delays where the coherence resonance is felt the strongest 
may be approximated by r = Xb/2 + n*Xb, with To being 
the period of coherent oscillations at r = 0. Noise- wise, 
with increasing r, D has to be adjusted to higher values 
to regulate the relaxation from the kink to the slow man- 
ifold while maintaining the entrainment to the proper 
frequency. In parallel, for stronger D, the representa- 
tion cloud of the firing fraction tends to disperse more, 
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FIG. 4. (a) and (b) show typical PPs of neurons par- 
ticipating the homogeneous global oscillations and clustered 
states, respectively. The latter are distinguished by a kink 
K, which is a signature of the D — r co-effect. The pa- 
rameter sets are r = 6, D = 0.0005, c = 0.1 in (a) and 
T = 2, D = 0.00025, c = 0.1 in (b). 

requiring a sufficient r for this effect to be averaged out. 

The interplay between D and r is further highlighted 
by exploring the behavior of the MF model (|2]). Local 
bifurcation analysis shows that the MF exhibits a suc- 
cession of super- and subthreshold Hopf bifurcations ^ , 
which account for the transition from the stochastically 
stable FP to the stable LC. Still, this scenario is con- 
fined to noise higher than here: analytical and numer- 
ical means corroborate the Hopf bifurcations to emerge 
about D « 0.0025 at relevant r. Now we argue that the 
approximate model is in qualitative terms able to cap- 
ture the clustering effect occurring for small D, c and r. 
Focus is on the finding that MF system predicts an on- 
set of cluster states by undergoing a global bifurcation 
for the parameter values around r = 2, D = 0.00025 and 
c = 0.08. At the given r and D, for c < 0.08 the approx- 
imate model has only the equilibrium, whereas around 
c c± 0.08 a large and a small LC are born via a fold-cycle 
scenario. Note how then the PP of the MF acquires the 
form qualitatively similar to those of the exact system's 
in Figs. [2fa) and[2jb). The two sections of the emerging 
MF orbit mimic the action of the fractions within the 
full population. This structure of the LC becomes unsta- 
ble under increasing c or r, i.e. for the stronger impact 
of the interaction term. Another interesting aspect to 
the approximate system is that it shows the complex LC 
to coexist with the FP, viz. the basins of attraction in 
Fig. E^b), which is a feature apparently absent in the ex- 
act model. However, the FP is located very close to the 
basins' boundary which indicates it to be stochastically 
unstable in the exact system for an arbitrary small noise. 

We have reported on a novel phenomenon where clus- 
tering within the homogeneous neural population is in- 
duced by an interplay of noise and time delay. This 
paradigm is distinct from most current explanations on 
how the clustered states may arise, for it does not treat D 
and r as destabilizing and detrimental, but rather as bi- 
ased toward the formation of dynamical structure in net- 
works that are unstructured both in terms of topology 
and local parameters. The analyzed model is minimal 
yet sufficient to display an interesting type of behavior, 



FIG. 5. (color online). Bistability in the MF model: (a) shows 
the trajectories converging either to the FP or the LC, de- 
pending on the initial conditions, whereas in (b) are displayed 
the two basins of attraction for r — 2, D — 0.00025, c = 0.1. 
possible only as an interplay of excitability, noise and in- 
teraction delay. Once the phenomenon is recognized as 
caused only by these qualitative properties one can study 
the effects of more realistic assumptions on the distribu- 
tion of neuronal properties and connection patterns. An 
interesting point concerns the derived MF model, which 
can aid in understanding the precise roles played by D 
and r. Notably, beneath the surface lies a more strati- 
fied phenomenon, where the subtle adjustment between 
the parameters affects the number of clusters, their con- 
figuration, stationary or dynamical character, as well as 
whether the cluster states occur monostable or coexist 
with the homogenous solution at the given population 
size. This framework could find application within the 

research on neural systems and other excitable media. 
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